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Shock formations are observed in granular avalanches when supercritical 
flow merges into a region of subcritical flow. In this paper we employ a 
shock-capturing numerical scheme for the one-dimensional Savage-Hutter 
theory of granular flow to describe this phenomenon. A Lagrangian moving 
mesh scheme applied to the non-conservative form of the equations repro¬ 
duces smooth solutions of these free boundary problems very well, but 
fails when shocks are formed. A non-oscillatory central difference (NOC) 
scheme with TVD limiter or WENO cell reconstruction for the conservative 
equations is therefore introduced. For the avalanche free boundary prob¬ 
lems it must be combined with a front-tracking method, developed here, to 
properly describe the margin evolution. It is found that this NOC scheme 
combined with the front-tracking module reproduces both the shock wave 
and the smooth solution accurately. A piecewise quadratic WENO recon¬ 
struction improves the smoothness of the solution near local extrema. The 
schemes are checked against exact solutions for (1) an upward moving shock 
wave, (2) the motion of a parabolic cap down an inclined plane and (3) the 
motion of a parabolic cap down a curved slope ending in a flat run-out 
region, where a shock is formed as the avalanche comes to a halt. 
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1. INTRODUCTION 

Snow avalanches, landslides, rock falls and debris flows are extremely dangerous 
and destructive natural phenomena of which the occurance has increased during the 
past few decades. Their human impact has become so significant that the United 
Nations declared 1990-2000 International Decade for Natural Disasters Reduction 
(IDNDR). Research on the protection of habitants from floods, debris flows and 
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avalanches is under way worldwide, and many institutions focus on the numerical 
prediction of such flows under ideal as well as realistic conditions. 

One of the models that has become popular in recent years is the Savage-Hutter 
(SH) avalanche theory for granular materials [33, 34]. In the past decade numeri¬ 
cal techniques were developed to solve the SH-governing differential equations for 
typical moving boundary value problems [33, 34, 13, 9, 10, 19, 16, 14, 8, 41, 6, 7]. 
These techniques are based on a Lagrangian moving mesh finite-difference scheme 
in which the granular material is divided into quadrilateral cells (2D) or trian¬ 
gular prisms with flat tops (3D). Exact similarity solutions of the SH-equations 
were constructed in spatially one-dimensional chute flows [33, 35, 31] and for two- 
dimensional unconfined flows [15, 12]. In the case of chute flows it was shown that 
the solutions obtained by the Lagrangian integration procedure approximate the 
exact parabolic similarity solution very accurately, and these theoretical and numer¬ 
ical results are in good agreement with experimental avalanche data. Similar agree¬ 
ment between theoretical, numerical and experimental data was also obtained for 
the two-dimensional flow configurations (cf. above references). In these Lagrangian 
schemes explicit artificial numerical diffusion was incorporated to maintain stabil¬ 
ity. In doing so the quality of resolution deteriorates. In fact, the adequacy of these 
numerical solutions can be challenged because of uncontrolled spreading due to this 
diffusion. It was also observed that the Lagrangian schemes loose their stability (or 
else unjustified artificial diffusion must be applied) whenever internal shocks are 
formed. This appears to occur whenever the avalanche moves from an extending 
to a contracting flow configuration. These shocks are travelling waves which form 
bumps with steep gradients on the free surface, which is thicker on the down slope 
side. It is therefore natural to develop conservative high-resolution shock-capturing 
numerical techniques that are able to resolve the steep surface gradients and iden¬ 
tify the shocks often observed in experiments but not captured by the Lagrangian 
finite difference scheme. 

The development of high-resolution shock-capturing schemes has a long history 
which we cannot even sketch here (see e.g. the classical references [3, 40, 11, 42] or 
the recent textbooks [25, 4, 20, 39]. The most common approach is to Hrst develop 
a one-dimensional TVD (total-variation-diminuishing) upwind scheme for a scalar 
conservation law and then apply it to systems using one-dimensional characteristic 
decompositions or approximate Riemann solvers. Upwind schemes have been used 
very successfully for gas dynamical calculations, where the Riemann problem can 
be solved exactly and many approximate Riemann solvers are available. For more 
complicated systems like the granular flow model considered here characteristic de¬ 
compositions are often not available, and the Riemann problem cannot be solved 
analytically. Therefore we have chosen an alternative approach to high-resolution 
shock-capturing, namely the recent non-oscillatory central (NOG) schemes first in¬ 
troduced by Nessyahu and Tadmor [30]. While upwind schemes are higher order 
extensions of the classical Godunov scheme, central schemes build upon the (also 
classical) Lax-Friedrichs scheme [23]. This scheme avoids characteristic decompo¬ 
sitions and Riemann solvers by the use of a staggered grid. When used together 
with piecewise constant spatial reconstructions, the Lax-Friedrichs scheme is more 
diffusive than Godunov’s scheme. However, when one combines the scheme with 
TVD-type piecewise linear reconstructions, it becomes competitive with the up- 
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wind schemes. Recently, central schemes have been extended in many directions, 
see eg. [1, 2, 18, 27] for multidimensional extensions, [32] for an adaptive staggered 
scheme, [28, 26] for third- and higher-order schemes and [22, 21] for central schemes 
on non-staggered grids, which are precisely at the borderline of central and upwind 
schemes. 

Here we adapt the second order NOC scheme of Nessyahu and Tadmor to include 
an earth pressure coefficient which has a jump discontinuity as the flow travels from 
an expanding into a contracting region, and to treat the source term which is due to 
the spatially varying topography and the gravitational force. The resulting scheme 
works well both in smooth regions and at shocks, which are captured within two 
mesh cells and without any oscillations. 

Besides the formation of shock fronts in the interior, avalanches may also have a 
vacuum front at their margins. Similarly as for the equations of gas dynamics, the 
hyperbolic system degenerates at the vacuum state. Many shock capturing upwind 
schemes produce negative heights at these points and subsequently break down or 
become completely unstable. While our NOC scheme is remarkably stable at the 
margins, it does not capture the vacuum front as well as the Lagrangian moving 
mesh scheme. To overcome this imperfection, we augment the NOC scheme with 
an algorithm that tracks the vacuum front. The combined front-tracking non- 
oscillatory central scheme is accurate and robust both at shocks and at the margins 
of the granular avalanche. 

The ensuing analysis commences in § 2 with the presentation of the governing 
SH-equations in conservative and non-conservative form; then the jump conditions 
of mass and momentum at singular surfaces will be stated and the solution to a sin¬ 
gle shock wave (a hydraulic jump) will be presented; § 2 closes with the construction 
of exact similarity solutions of a parabolic heap moving down a rough incline. § 3 
introduces the numerical techniques; at first the Lagrangian integration technique 
is described; it is followed by the presentation of the non-oscillatory central (NOC) 
scheme. In §4 we augment the NOC scheme (which uses a fixed Eulerian grid) 
with a Lagrangian type front-tracking method in the marginal cells. § 5 elaborates 
on numerical results. The travelling shock wave cannot be handled by the La¬ 
grangian method, but the NOC scheme can do so with very little diffusion accross 
the shock. On the other hand, the parabolic similarity solution is well produced 
by the Lagrangian integration technique, but much less accurately by the NOC 
schemes unless the Lagrangian front-tracking is introduced for the marginal cells. 
It is also shown that the NOC scheme with piecewise linear spatial reconstructions 
applying standard TVD type slope limiters exhibits some oscillations near smooth 
local maxima. We remove these oscillations by incorporating a piecewise quadratic 
WENO (weighted essentially non-oscillatory) reconstruction into our scheme. Our 
final numerical experiment combines all the difficulties treated in the paper: an 
avalanche with a vacuum front at the margins expands as it flows downhill and 
contracts as it hits the flat runout (so the earth pressure coefficient changes discon- 
tinuously at the transition point). As the avalanche comes to a halt at the bottom, 
a shock wave develops and propagates upslope. Our NOC front tracking scheme 
handles this challenging flow very satisfactorily. § 6 presents conclusions and gives 
an outlook to further work. 
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2. GOVERNING EQUATIONS 

A detailed derivation of the Savage-Hutter theory has been given in [33, 34]. Here 
we confine ourselves to a brief description. Although cohesionless granular materials 
exhibit dilatancy effects numerous experiments have confirmed that during rapid 
dense flow it is reasonable to assume that the avalanche is incompressible with 
constant uniform density po-” During flow the body behaves as a Mohr-Coulomb 
plastic material at yield. As the avalanche slides over the rigid basal topography 
a Coulomb dry friction force resists the motion. The basal shear stress is there¬ 
fore equal to the normal basal pressure multiplied by a coefficient of friction tan <5, 
where 6 is termed the basal friction angle [19]. Scaling analysis isolates the phys¬ 
ically significant terms in the governing equations and identifies those terms that 
can be neglected. Plane flow configurations are our focus in this paper, so depth 
integration reduces the theory to one spatial dimension. The leading order, dimen¬ 
sionless, depth integrated equations for the local thickness of the avalanche h and 
the momentum hu {u is the downslope velocity) reduce to 

dh 9 

_ + _(ft„)=0, (1) 

+ = ftSi- (2) 

with net driving force 


Sx = sin^ — sgn(M) tan(5(cosC -I- Xku^) — scos(-^, (3) 

where x is the arc length measured along the avalanche track, denotes the height 
of the basal topography relative to the track (usually = 0 in one spatial di¬ 
mension) and C and Xk are the local slope inclination angle and curvature of the 
track, respectively. The term sgn(t6) selects the orientation of the dry Coulomb 
drag friction, and £ <C 1 is the aspect ratio of a typical thickness and length of the 
avalanche. Note that equations (1) and (2) are written in conservative form [8], 
while in the original SH theory the smoothness assumption allows the momentum 
balance equation to transform to an evolution equation for the velocity, viz., 


dt6 

dt 


dh _ I dj5x 
' dx 2 dx ' 


( 4 ) 


The factor /3x is defined as (5x = e cos(Kx and the earth pressure coefficient Kx is 
given by the ad hoc assumption 


Kx = 


^Xact 

for 

dujdx > 0, 

^Xpass 

for 

dujdx < 0, 


( 5 ) 


with 

= 2 (^1 T \/l-cos2 (/)/ cos2 6^ sec^ </> - 1, (6) 

and (j) is the internal friction angle of the granular material. Note that the values of 
the earth pressure coefficient Kx are based on the postulation of a Mohr-Coulomb 
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plastic behaviour for the cohessionless yield on the basal sliding surface, see Savage 
& Hutter [33, 34] for details. In this theory the earth pressure coefficient Kx is 
assumed to be function of the velocity gradient, i.e. Kx = Kx{du/dx). 

The governing equations look like the shallow-water equations, but because of the 
jump in the earth pressure coefficients the source term Sx and the free 

boundary at the front and rear margins, it becomes much more complicated to de¬ 
velop an appropriate numerical scheme to describe the flow. The original Lagrange 
finite-difference scheme [33] is implemented for the equation system (1) and (4) in 
Lagrangian form, with primitive variables h and u. The shock capturing scheme 
developed here is applied to the system in conservative form (1) and (2), where 
the conserved quantities are the avalanche thickness h and the depth integrated 
momentum m = hu. 

In vector notation, equations (I) and (2) take the form 

wt + fx= ( 7 ) 


where 


’"=(„)• I = (mVft+“ ‘={hs.)- W 

This form is more convenient for mathematical analysis than (1) and (2). 

2.1. Jump Condition and Travelling wave 

The Savage-Hutter theory can be used to model the upslope propagating travel¬ 
ling shock wave observed in experiments [5, 7] by introducing the jump conditions 
(see Fig. 1) of the balance equations (1) and (2) for mass and momentum 


|/i(M-y„)i = o, (9) 

{hu{u - Vn) + iPxh"^} = 0 , (10) 

where 14 is the normal speed of the singular surface. Let us suppose that [(3x\ = 0 

(for example, this is always satisfied \i (j) = 5, i.e., = ^xpass)- Substituting 

(9) into (10) (i.e. eliminating 14) yields the following relation between the depth 
ratio, H := h~/h'^, and the velocity difference 

. (11) 


For an upslope travelling shock wave with travelling wave speed 14 and corre¬ 
sponding depth ratio H, the factor /Sx is a function of material and topographic 
parameters, </>, S and 4 which are given by the selected material and topography. 
Provided that the depths before and after the shock, and h~, are known (they 
can be determined by experiment) and the downslope velocity is also given (it is 
normally equal to zero), then the upslope velocity can be determined by using (11) 


= u ± 


H- 1 


H 


Pxh' 


. 77+1 


-|l/2 


( 12 ) 
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h+ 


a 


-Vn 




FIG. 1. The plane travelling shock wave can be interpreted as a jump in thickness and 
velocity separating the body of the avalanche into two parts on a plane with inclined angle 
and h~ are the thicknesses of both sides and u"*" and u~ are the velocities, respectively, whereas 
this jump travels with velocity Vn up-slopes. 


Note that the term under the square root is positive for all positive H. If H = 1 
then u~^ = u~ ^ which indicates that no shock wave (discontinuity) takes place. 
Thus, velocity jumps and depth jumps occur together. 

By inspection of the mass balance equation (9), the velocity of the shock is given 

by 




Hu — u'^ 
H -1 


u T 



H + 1 


1/2 


2i?2 


Note that as tends to h = h , then 11 + tends to u = u and 


(13) 


14 u =F [iSxh]^^'^ , 


(14) 


so we have recovered the characteristic speeds of the shallow water equations. Now 
we apply Lax’ shock inequalities [24] to single out the physically relevant branches 
of the shock curves: for the first family, with characteristic speed u — y /we 
require that 


u+ - [(3xh+] > Vn = u- 



H + 1 


1/2 


2iL2 


> u 




which implies H > 1 (recall that the upslope state “+” lies to the left of the shock). 
Analogously, for the second family, with characteristic speed u + ^f3xh, we obtain 
H < 1. For example, an upward jump {h~^ < h ) can only be carried by a shock of 
the first family, and in this case u~^ > u~ > Vn, so particles which cross the shock 
are condensed and slow down. 


2.2. Similarity Solution 

Consider the motion of a finite mass of granular material along a flat plane, i.e. C 
is constant and Ak = 0 in (3). In [33] one particular similarity solution to a moving 
boundary problem of finite mass was derived; this solution is now generalised (see 
[36]. To this end we introduce a moving coordinate system with velocity 

Mo(i) = Mo(0) + / (sin C — tan (5 cosC) dt 

Jo 


(15) 
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on an plane with inclination angle This velocity is due to the net driving force 
Sx in (3), where we assume that the velocity is positive for positive times, i.e. 
sgn('u) = 1. The relative velocity u in the moving coordinate system is then given 

by 

u = u — uo{t) . (16) 

A symmetric bulk is considered and the origin of the moving coordinate system is 
selected to lie at the centre where the surface gradient, dh/dx, is zero. To keep the 
symmetric depth profile during the motion the relative velocity is further assumed 
to be skew-symmetric, u(^, t) = t), where 


^ = x- [ Uo(t')df (17) 

Jo 

indicates the distance from the origin in the moving coordinates. Provided that 
g{t) is the distance from the coordinate origin to the margin at time t, the physical 
domain occupied by the granular mass can be mapped from [—g{t), g{t)] to the 
fixed domain [—1, 1] by 

?7 =y uo(f')dt'| , where 7 ?g[-1,1]. (18) 

With this coordinate mapping, {x, t) —>■ (r/, r), the model equations (1) and (2) 
reduce to 


dh 9'nh la 


(19) 


du g' du 
dt ^ g dg 


I f^du a 


= 0 , 


( 20 ) 


where the r is again replaced by t and we have used g' = dg/dt = —uo/g. 

Now we assume that u{g, t) varies linearly in g. Since the margins move with 
relative speeds ±g'{t), this yields u{g,t) = gg'{t). Now the evolution equations (19) 
and (20) reduce to 


9 


99 


dt ' g 

, 

9 dg 


h = 0, 

= 0 , 


( 21 ) 

( 22 ) 


where g" = d?g/dt^. Integrating (22) subject to the boundary condition either 
h(g = 1) = 0 or h{g = —1) = 0, it follows that the thickness is described by 


h{g, t) 


9{t)9"it) 

2/3. 


(l-g^). 


(23) 


This implies that the avalanche body keeps a parabolic thickness distribution during 
the motion. With the thickness distribution (23) one can easily obtain the total 
mass M to be 


M = 



t)d^ 


Hv, t)g{t) dg 


2 9 " 9 ^ 

3 Px 


-1 


(24) 
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Since mass is conserved, 


(25) 


This relation can also be derived directly from the mass balance equation (21). 

Changing the independent variable t to g{t) and letting p{t) = g'{t), equation 
(24) can be written as 


dp K 

dp g^ 


(26) 


where 3f3xM = 2K. The similarity solution is then obtained by solving (26) with 
initial condition, p(0) = go and p(0) = po 


p‘^{t) = 2K 




(27) 


With the definition ctg = Pg = Po ^^d G = (og + /3g )g it follows that 


VGG' 

y/G-2K 


(Og + 


(28) 


We now use the relation 


d 


VG^/G - 2K + 2K \a{VG + y'G - 2K) 


Vg 

y'G-2K 


and integrate equation (28) to yield 


VGVG - 2K + 2K ln{VG + VG - 2K) 


VGy'G - 2K + 2K \n{VG + VG - 2K) 


t=0 


{ag+Pgf/'^t. 


With Po = 1, Po = 0 we obtain the Savage-Hutter solution [33] 
Vg^/g- 1 + In (Vs + \/p - l) = C 


(29) 


(30) 


for which g{t) > 1. Both (29) and (30) are implicit evolution equations for g{t). 
Once g{t) is deduced, with the presumption u{p,t) = pg'{t), the complete solution 
is then given by (23) and (27), 


u{r], t) = p 




Hv, t) 


3M 

4^ 


(1-p^), 


(31) 


where g is defined in (18). In the present similarity solution it is presumed that 
m/|u| = 1, which means that u > 0 for all t > 0. From (16) and the presumption 
u{r],t) = r]g'{t) it follows that 


u{t) = uo{t) + u(t) >0 g'{t) < uo(t), for all t > 0. (32) 


It is very important to verify that the velocity is consistent with condition (32) 
to keep the parabolic similarity solution valid. The generalisation (29) of (30) 
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was needed to have exact solutions with non-vanishing initial velocities (for further 
details see [36]). 


3. NUMERICAL SCHEME 

The numerical schemes employed in this paper are designed to explicitly solve 
the system of equations in ID and we here introduce a Lagrangian algorithm and 
an Eulerian shock-capturing NOC (Non-Oscillatory Central) scheme. 

In the Lagrangian technique [33, 34] the avalanche body is divided into several 
cells. The purpose is to find the velocity of the cell boundaries in order to deter¬ 
mine the cell boundary locations for each time step; so it is a moving-grid method, 
whereas, the NOC scheme is built on a stationary uniform grid and gives a high 
resolution of the shock solutions without any spurious oscillations near a disconti¬ 
nuity. 

In the Lagrangian method the value of the depth /i" is defined as the volume 
average within the cell for time t", which is bounded by &j_i(t) and bj{t), and 
the boundary bjit) moves with the velocity uj. Whilst, in the NOC scheme the 
value of the discretised variable C/", U = h, m is defined on the mesh as the volume 
average within the mesh cell centred at position Xj for time where the 
cell is bounded by Xj+if 2 and Xj-if 2 - 


3.1. Lagrangian method 

In the Lagrangian method [33, 34] the avalanche body is divided into N material 
cells, where x = bj-i{t) and x = bj(t) denote the boundaries of the cell j at time t, 
see Fig. 2. These boundaries move with the avalanche velocity, i.e. 

=Uj{t) =u{bj{t),t). 

Integrating the mass balance equation (I) over the cell yields 


/ 


C-i 





hdx = 0 



(33) 


and implies that the volume (mass) of the cell is conserved during the motion. 
Because of this, the mean height of the cell can be determined by 


j^n _ ^cellj 

^ “ 6 ” - ' 


(34) 


The computations proceed as follows. It is assumed that 6”, /i” and 

•J _ 1 _ 1 ^ 

are given as initial values and the new location of the cell boundary b^ after an 
elapsed time At is given by 




+ At ti 


n+l/2 

3 


(35) 
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FIG. 2. The avalanche body is divided into N elements with average depth hj, where cj 
is the centre of the element. 


Note that here the velocity Uj indicates the boundary velocity of bj . The momentum 
balance (4) allows the velocity of the cell boundary at time to be determined, 


„+l/2 n-l /2 , n A f f dh'^ ^1+1/2 f d{cOsCK^)\‘\ 

The net driving acceleration s" as given by (3) is 

sin^j — sgn^u””^^^^ tan5 |cosCt + Xkj i ~ ^ 


= 


dz'’ 

dx 


(36) 

(37) 


where ((,■ represents the local inclination angle, Kj is the local curvature, and 
denotes the local basal topography. Note that the last term at the right-hand side 
of (36) contains the gradient of the earth pressure coefficient, which is neglected in 
the numerical scheme of Savage and Hutter [33, 34]. 

The earth pressure coefficient Kx is determined by the ad hoc definition 


{KxT, = 



for 

n—ll2 

^+1 

> 

n-l/2 

Kx 

'^pass * 

for 

n—1/2 
^j+i 

< 

n-l/2 

^3 


(38) 


in [33, 34]. The surface (depth) gradients in (36) are determined by the depths of 
the adjacent elements 


\dx)^ ^ ^ 

where c” represents the centre of the cell, c" = (5” -|- b^_-^)j2, at time t = t", 
see Fig. 2. The height at the cell boundary, /ij+ 1 / 2 , is given by their mean values 
in adjacent cells, /ij+ 1/2 = ^ [hj + hj+i ), and the gradient of the earth pressure 
coefficient is 


/ d{cos C,Kx) Y 

V dx 


(40) 


However, while this method is excellent for classical smooth solutions, it looses 
numerical stability if shocks develop. Shocks are initiated when the avalanche 
velocity is faster than its characteristic speed and the avalanche front reaches the 
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a) 


b) 


J+i 


J+2 



j+112 

»-< 

j+312 




- •—^ 


^ X n 


j+112 


j+114 


J+314 


J +1 


A2 


j+1 


FIG. 3. Diagram of NOGS scheme, a) Grid points computed in the NOGS method, 
b) NOGS computational diagram, where • indicate the grid points at time level n and n -|- 1, 
■ represent the quadrature points for the fluxes f across the cell boundaries, <> the quadrature 
points for the source terms s and A those for the staggered cell averages at the original time i". 


base of the slope or a solid wall. Many detailed investigations about granular shocks 
were made by Gray and Hutter [5], in which the shock waves are concerned to be an 
important property in the granular flows. To avoid the numerical instability caused 
by the shocks, an artificial viscosity term ^d'^ujdx^ is introduced and added to the 
right hand side of (37) for numerical stability, e.g. [33, 34] and [14], where the 
artificial viscosity /i was found to have values between 0.01 and 0.03. 

3.2. NOC Scheme 

The Non-Oscillatory Central Differencing (NOC) scheme of Nessyahu and Tad- 
mor [30] is a second order accurate extension of the classical Lax-Friedrichs scheme 
[23]. Let us briefly review the NOC scheme: 

We consider the Savage-Hutter equations in the conservative form (7), (8) with 
w = {h, nif" as basic variables. Let w” denote the cell average over interval 
at time t", and let 


w(x,r) = w7 + ^^w' (41) 

be a piecewise linear reconstruction over the cell, where W' denotes the cell mean 
derivative determined by a TVD-limiter [25] or a central WENO cell reconstruction 
[26]. The main conceptional difference between the NOC schemes and standard 
upwind finite difference schemes is the use of a staggered grid. At time = 
t”-|-At, the cell averages w"”])! are evaluated over the intervals [xj,Xj+i], see Figure 
3. As a consequence, the boundaries of the cells at the new time level are the centers 
of the cells at the old time level, namely the points Xj and Xj+i. At these points, 
the piecewise polynomial reconstruction (41) of the cell averages at the old time 
level t" is smooth, and it remains so for t < under an appropriate restriction 
of the timestep (see (49) below). Therefore, the flux across the boundaries of the 
cells at the new time level may be evaluated by Taylor extrapolations using the 
differential equation and standard quadrature rules. Here we use the midpoint rule 
in time to achieve second order accuracy. The resulting update takes the form 


3 + 1/2 




j + l/4 


C+3/4 


Ai f ^^+1/2 ^71-1-1/2 

a7 \^N-i-i ~ N 


At / 
2 


,+ 1/2 

’i+1/4 


'i+l/2\ 
"i+3/4 ) 


(42) 













12 


Y.C. TAI, S. NOELLE, J.M.N.T. GRAY & K. HUTTER 


as illustrated in Fig. 3^. The values of and are determined by the 

reconstruction (41) over the and (j + 1)*^ cell, i.e. 

= wj + iw', w"+3/4 = w”+i - iw'+1. (43) 

The transport flux f at the quadrature points {xj, i"+i/2) and (a^^+i, is 

approximated by Taylor extrapolation in time, 

f„+i/2 ^ f ^ (aw/at)”, (44) 

and similarly, the source terms s at the quadrature points {xj+ 1 / 4 , t"+i/2) and 
(a;j_|_ 3 / 4 , are approximated by space-time Taylor extrapolation 


n-i- 1/2 

Sj+1/4 = s I w 


0 + 1/4 ) 


_n+1/2 _ n 

'^j+ 1/4 - 


+ ^ (aw/at)” -b tw'. 


n+1/2 

®J+3/4 


= s [w 




i+l/2\ 

'j+3/4 ) 


(45) 


—n+1/2 —I At / o—/ 1—/ 

w^-+3;4 = w^.+4 + — (aw/at)^.+4 - 


The temporal derivative (aw/at)" in (44) and (45) is determined by using (7), 

{dw/dtJJ = - (af/aa;)" +s” = -A^ w'/Aa;-b s”, (46) 

where 


(df/dx); = (A); (aw/acr);, a = af/aw = ( L ) (47) 

and A is the Jacobian of f. Alternatively, one may also use the Jacobian-free 
approach of Nessyahu and Tadmor [30] and set 


(af/aa:); = f'/Ax, 

where the cell mean derivative f' of the flux is again determined by a TVD-limiter. 
Let be the maximum wave speed, 

^max _ « / , Uj = rrij/hj for hj 0 . (48) 

all j \ 

The CFL condition 

|+a-“"|<i, for all j (49) 

is needed to guarantee that the solution remains smooth at the space-time quadra¬ 
ture points, so that the Taylor expansions (44) and (45) are justified. 

Note that the NOC scheme (41) - (49) completely avoids the expensive Riemann 
solvers used in standard upwind schemes on non-staggered grids. The resulting 
staggered schemes are easy to code, computationally efficient and can be applied to 
general systems of conservation laws, where the solution of the Riemann problem 
(i.e. the initial value problem with piecewise constant data) may complicated or 
even impossible. 
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4. FRONT-TRACKING METHOD 

In many applications, the region covered by the granular material has a finite 
extension and is limited by a free boundary which moves with the flow velocity. 
Outside this region, there is vacuum, so the avalanche height h and momentum m 
are zero, and the velocity u = m/h is not well-defined. The Lagrangian method 
handles this situation automatically, since the computational domain moves with 
the material flow. The NOG scheme discretizes the differential equations on a 
stationary uniform mesh. Note that in general the margin points (the front 
margin) and (the tail margin) lie between grid points, so that it is impossible 
to point out the margin locations without extra treatment. Furthermore, it is not 
straightforward to determine the proper cell reconstructions over the margin cells. 
Fig. 4 illustrates an example of depth reconstructions over the front margin cell 
determined by various TVD limiters. Here and in the following we suppose that 
at time the front margin lies in the cell, and the tail 

margin in the cell, Xi_i < x^; < X(_|_i. 





FIG. 4. Example of the depth reconstruction (solid line) determined by different TVD 

limiters, where the circles denote the cell average. The front margin lies in the cell. In the 
Eulerian scheme one cannot determine where the margin lies. Outside the margin there is no 
material, so that the average depths of the cells / + ^ > 1 are equal to zero. Different limiters 

lead to different outflows from the avalanche body. 


Since our quadrature rule for the fluxes (44), (46), (47) uses a Taylor expansion 
of the solution, different limiters will lead to different values of the integrals of the 
fluxes across x/ x and Xt x To complicate the situation even 

further, part of these boundaries may lie in the vacuum region. Note that the fluxes 
across these boundaries determine the outflow from the avalanche body, so non- 
appropriate cell reconstructions over the margin cell may cause too much outflow 
from the avalanche body or even result in a negative depth around the margin, 
see Fig. 5o. Thus, the difficulty is not only to determine the correct numerical 
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flux at the grid point Xf, the wrong numerical flux may also cause vast stability 
problems. Adding a thin layer over the whole computational domain can circumvent 
the numerical stability problem, but it is then difficult to determine the locations of 
the margins, and the numerical flux out of the avalanche body would even become 
unexpectedly large, which results in large numerical diffusion, while there will be 
permanent outflow from the avalanche body. Therefore, a more refined treatment 
is needed for the evolution of the avalanche margins. 

In [29], Munz developed a method to track vacuum fronts in gas dynamics. His 
approach is based on appropriate reconstructions of cell averages behind the front, 
and the solution of a vacuum Riemann problem which is used to track the mar¬ 
gin locations at every time step. Here we develop an alternative front tracking 
method, which is based on a piecewise linear spatial reconstruction of the conser¬ 
vative variables up to the front and Taylor extrapolations in time. Contrary to [29] 
our approach is Riemann-solver free and therefore fits perfectly into the framework 
of central schemes. 

The structure of our front-tracking algorithm is as follows: At the beginning of 
each timestep (at time t„), the cell averages w” of the conservative variables and 
the position of the margin points x'^^ (front) and xjj (tail) are given. In the first 
step, a piecewise linear reconstruction of the data is defined, the front (tail) velocity 
is determined and the front (tail) is propagated from time to t^+i- In the second 
step, the conservative variables are updated via 


n+1 

i-i/2 


1 

Ax 


w(a;, t") Ax 


'Xj-i 



i{xj-i,t)} At 


1 

+-7 — / / s{x,t)AxAt, 

Jtn 


(50) 


Away from the front, the integrals are evaluated by the midpoint rule as in (42). 
Special care has to be taken in the two margin cells (the cells containing the front 
and the tail). Each of the integrals on the RHS of (50) may contain parts of the 
vacuum region. Therefore, we need to replace the midpoint rule by more delicate 
quadrature rules over the region covered by the granular material. 

In order to guide the reader through the details of the algorithm, we would like 
to give an outline of the rest of this section. In Section 4.1, a particular piecewise 
linear reconstruction of the conservative variables near the front is derived. In 
Section 4.2, the front velocities are computed, and the fronts are propagated to 
the new time level. In Section 4.3, four cases are distinguished for the location of 
the front relative to the fixed underlying grid, and their geometry is discussed. In 
Sections 4.4, 4.5 and 4.6, the three integrals on the RHS of (50) are treated: the 
data, the fluxes and the source terms. In Section 4.7, a special space-time Taylor 
extrapolation of the conservative variables near the front is derived, which is needed 
to compute the solution at the space-time quadrature points of the three integrals. 
Section 4.8 summarizes the algorithm. 


4.1. Reconstructing the conservative variables 
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Suppose as before that the front margin is contained in the cell, and the rear 
margin in the cell, 

Xpt G 

^Ti G 

We require that the piecewise linear reconstruction w(x,t„) satisfies the following 
two criteria: 


• first, it should vanish at the margin points, and 

• second, it should preserve the cell averages. 

These criteria uniquely determine the reconstruction in the margin cells. If we 
denote the cell averaged depths of the front margin and rear margin cells hy hf and 
ht, the depth reconstruction is defined by 


hf{x)=a’f{x-XFt); cr’f = 

for the front margin cell and by 
htix) = at(x -XTi); <Jt = 


—2 hf Ax 

{xFt - Xf_il'2?' 


2 ht Ax 

{Xt+l/2 - XtiY' 


for X G [a;/_i/2, xpt], 


for X G [xTu Xt+1/2], 


(51) 


(52) 


for the rear margin cell. Outside the margin the depth is equal to zero. The xpt 
and Xf-i /2 represent the locations of the front and the internal boundary of the 
front margin cell, respectively (see Fig. 5). The xti and a;t+i /2 denote the locations 
of the rear and the internal boundary of the rear margin cell, respectively. The 
reconstruction of to = hu, m{x), is defined analogously. 


4.2. Propagating the front 

Besides being very natural and simple, our definition of the reconstructions over 
the margin cells has the advantage of leading to uniquely defined, constant values 




FIG. 5. The reconstruction of the depth hf(x) within the margin cell, (a) Cell 

reconstructions based on TVD-limiters cannot determine the location of the margin point. Non- 
appropriate reconstructions over the margin cell may result in wrong values of the flux at the 
gridpoint Xf, which may cause too much outflow from the avalanche body, (b) Our front tracking 
method uses the unique piecewise linear reconstruction hf(x) over the margin cell which vanishes 
at the margin point xpt and preserves the cell average. Thus, a reasonable flux at xy is expected. 
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4.3. Intersecting the front and the grid 

Once the new location of the margin is given, the new margin cell at the next time 
step is then determined. The CFL condition (49) guarantees that At < 

Ax/2, so the margin point xptjTi can at most pass through gridpoint Xf/t during 
one time step. For example, with this condition the front can only lie in either 
one of the two adjacent cells of the margin cell, which are the (/ — and 
y(/ + cells, see Fig. 6. There are four possible cases for the motion of the front 
margin point: 

• case I : < a;/, and x/_i < x'^^ ^ a;/, 

• case II : Xf < x'^^ < X/+1, 

• case III : 

^ — ^/+1/2 5 

• case IV : x^( > Xf, and xy_i/2 < x^Y^ ^ a;/, 

where x^^ and x^^^ are the front locations at t” and respectively. In cases I 
and II, the front does not pass gridpoint x/, while in the cases III and IV it does, 
see Figure 6. In each case we have to determine the cell averages of the relevant 
cells and by integrating the governing equations over [x/_i,x/] x 

and [x/,x/+i] x respectively, i.e. we have to evaluate the three 

integrals on the RHS of (50). These integrals involve the data w, the fluxes f and 
the source term s. In the following, we derive quadrature rules which are exact for 
linear functions. The tail margin can be treated completely analogously. 


4.4. The integral of the data 
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First we integrate the linear reconstruction w(a;,t„) of the data at time over 
the interval [xf-i,xf]. In cases I and III, this interval contains the front, while it 
does not in cases II and IV. We obtain 


w(.,r)dx = 


5 Wj_ 3/4 + Wj in cases I and III 
5 (wj_ 3/4 + in cases II and IV 


(56) 


Here Wf_ 3/4 is given by (43), and by (51). With a given front location 

and Wj it is 

21 


= 1 - 


^Ft 

^Ft - ^f-1/2 


(57) 


Next we consider the integral over the interval [x/,x/+i]. Using (51) once more we 
obtain 


1 

Ax 


w(x, t") dx 


/Xf + i 


0 in cases I and III 


yx'^^—Xf-i/2 J 


in cases II and IV 


(58) 


4.5. The integral of the fluxes 

Due to the restriction of the timestep, the only grid-position which is possibly 
intersected by the front during the time-interval [t„, tn+i] is x = x/. Therefore, the 
flux at Xf-i can be evaluated exactly as in the interior of the domain, 

^ f(w(x/_i, t)) dt = (59) 


where is given by (44). The flux at x/+i vanishes, since this point lies in 

the vacuum region during the whole time interval. It remains to compute the flux 
at X/. In cases III and IV, where the front crosses x/, we use the midpoint-rule in 
time over that part of the interface which lies within the region covered by granular 
material. Let i and At be the midpoint and the length of this time interval. If t* 
is the time at which the front intersects x/, defined by 


+ {t* ~ tn)Upf — Xf, 


then 


and 


t = 


{tn+i +t*)/2 in case III 
{tn +t*)/2 in case IV 


At = 


tn+l -1* 

[ t*-A 

The midpoint rule for the flux now gives 


in case III 
in case IV 


1 

At 


f(w(x/, t) dt 

.+1 


0 in case I 
‘ m case 11 

in cases III and IV 


(60) 


(61) 


(62) 


(63) 
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Here = f(w(a;/,f)). In Section 4.7 below we will extrapolate the solution w to 
the quadrature point {xf,t). 


4.6. The integral of the source term 

The source term s has to be integrated over the quadrilateral regions shown in 
Figure 6. Let us call these areas of integration H. In the following lemma, we give 
a quadrature rule which is exact for linear functions vanishing at the front. 

Lemma 4.1. Let a,b,T > 0 and 

Ll := {{x,t) : i < t < i + T,x < X < X + a + {b — a) {t — i)/T}. 


Let s be a linear function over LI which vanishes at the boundary x = x + a+(b — 
a){t — i)/T. Then 


J j s{x,t)dxdt 


1 a‘^ + ab + b^ 
3^ a + b 


s(x, i + t /2) 


ujs{x, i + t/2). 


(64) 


Proof: W.l.o.g. let x = t = 0. The general form of s is given by 
s(a;, t) = {x — a — {b — a)t/T)a 

where ct is a real constant. W.l.o.g. let cr = 1. Then a direct computation gives 
that 

J J s{x,t)dxdt =—^{a^ + ab + b'^) = ujs{0,^) 

n 

□ 

Equation (64) may be interpreted as a special quadrature rule with node {x,t + 
t/2). We have choosen this node because it appears also in the quadrature rule for 
the fluxes treated in Section 4.5 above, so we can minimize the evaluations of the 
solution w. 

In the following we apply the lemma to the four cases. Let LI be the region covered 
by the granular material. First, we compute the integral over the intersection of LI 
with the union of the (/ — 1/2)*^ and the (/ -I- 1/2)*^ cell, H fl ([a;/_i, x/-|-i] x 
[t”,f”+^]). Using X = Xf-i, i = tn, a = — x/_i, b = Xp/i^ — x/_i and t = At 

in Lemma 4.1 gives 

(65) 

with 


UJf-i = 


n ^f+i 

s(x,t)dxdt = 

7-1 

At {xfrt - Xf-i)"^ + {x^t - Xf-i){Xp+^ - Xf-i) + {Xp+^ - X/_i)2 


+ ^Ft 2x/_i 


( 66 ) 


Similarly, for the integral over H fl ([x/,x/+i] x [t”,t”+^]) we obtain 

s{x,t) dxdt = w/Sy, 


(67) 


'Xf 
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where t = in cases I and II and t is given by (61) in cases III and IV, and 

the wheight is given by 


0 


in case I 




At i^pt ^f) Xf) 




-2xf 


in case II 


— in case III 

^{xp^ — Xf) in case IV 


( 68 ) 


Here At is given by (62). The integral over [xf-i,Xf] x is then computed 

by substracting the integral over [xf,Xf+i] x [t”, t”+^] from that over [xf-i,Xf+i] x 


I I s(a;, t) dxdt = — w/Sy. 

This completes the definition of the quadrature rules for the three integrals on the 
RHS of (50). It remains to extrapolate the solution w to the new quadrature point 
ixf,t) near the front. 


(69) 


4.7. Determination of the physical quantities at t 

In cases III and IV the margin point passes the cell boundary Xf at t* and goes 
into the neighbouring cell. The outflow in case III and the inflow in case IV through 
the cell boundary at Xf as well as the source term in the new and old margin cells 
are essential for determining the cell average of the margin cells in the front-tracking 
method. 

In case III the physical quantities flow through the boundary Xf into the {f + ^Y^ 
cell during the time interval [t*, The outflow is approximated by the value 

at {xf,t), where t = i(t"+^ + t*), which is determined by using Taylor series 
expansion at the point (xp^, t„) with respect to space and time. Using the margin¬ 
cell-reconstructed slopes and the mass balance equation (I), the avalanche depth h 
at {xf, t) is then given by 


f ~ + dx 


(Xpitn) 


(xf - x%) + ^ 


(t-r) 


dx 




{xf -x'Yt)- 


dm 

dx 


ii-n 


(Xpitn) 


(70) 


= (aTfixf-x-pY-i<xn]ii-e): 


= {a^rf [{xf-x-pY-u-pYi-n], 

where (cr^)^ and (cr™)^ are the slopes of the margin cell reconstructions defined in 
§4.1, and we have used the relation (cr™)^ = w^((cr^)y. Similarly, the momentum 
is approximated by 


f 71 


(x’^ttn) 


(xf -xYt) + 


dm 

'W 




{i-n. 


(71) 
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Using the cell reconstruction of the margin cell and the momentum balance 
equation (2), equation (71) becomes 


dm 




{xf - - 


5 (^ + 4 /?.) 


dx 


(t-r) 

(^Ft 


/ 2m 9m m^ 9/i 


/3a;h 


(9h ddx 


^ h2 9x ' "'"'"9x ' 2 9x 

= {a-^Tfixf - xl,) - {i - t") 


(t-r) 


{x'p^itn 


= Upt h). 

(72) 

In case IV the physical quantities at the boundary (x/, t) are determined in the 
same way, but the time points are defined differently: t = (t"+^ -|-t*)/2 for case III 
and t = + t*)/2 for case IV. 


4.8. Summary of the front-tracking algorithm 

The front-tracking algorithm may be summarized as follows: 


Here 


T77^+l 

W/-1/2 


At 


At 


- + (1 - a/)w/ - -^{f + —f/- 


^+1/2 


Ax 


W /_1 „+i /2 


Ax ^ 


'^/+l/2 


afMVf 


Ax ^ Ax ^ 


Uf 


At 


t 


0 in cases I and III 

n \ 2 

— I in cases II and IV 

Ft ^/-l/2/ 

0 in case I 
At in case II 
tn+i — t* in case III 
. t* — tn in case IV 

tn in case I 

t„+i /2 in case II 
(t„+i -I- t*)/2 in case III 
. (tn + t*)/2 in case IV 



(73) 

(74) 

(75) 

(76) 

(77) 

(78) 


The wheights w/_i and w/ are defined in (66) and (68). The values of w(x/,t), 
needed to determine f| and s^, are defined in (70) and (72). This completes the 
definition of the update at the front margin. The tail can be treated completely 
analogously. 
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5. NUMERICAL RESULTS 
5.1. Travelling shock wave 

In this test problem we are concerned with granular flow on a plane {Xk = 0) 
inclined chute (0 < a: < 36 dimensionless units), where the internal and basal 
friction angles are both presumed to be equal to the inclination angle, (f> = S = 
( = 40°. That implies a non-accelerative flow, Sx = 0, whose earth pressure 
coefficient is constant Kx = Kx^^t = Selecting e = 1 and using (6) yields 

/3x = ecosl^Kx = 1.84477. A jump of thickness H = h~/h'^ = 3 with h'^ = 
0.3, h~ = 0.9 is presumed at a; = 24. By virtue of (12) the velocity difference is 
then determined, u~^ — u~ = 1.2148317, where the positive sign is selected. Since 
an instability was expected close to m = 0 as a singularity by sgn(u), the downslope 
velocity is assumed to be u~ = 0.1, so that the term sgn(M) is always unity. The 
initial condition of this test problem is defined as follows: 


/i(x, 0) = j 

f 0.3, for 0 < X 
( 0.9, for 24 < 

< 24, 

X < 36, 

(79) 

■u(x, 0) = j 

f 1.3148317, for 
( 0.1, for 

0 < X < 24, 

24 < X < 36. 

(80) 


From (13) the velocity of the upslope travelling wave is then expected as 14 = 
—0.50741585. For the boundary condition a constant inflow at x = 0 and an 
outflow condition for x = 36 are introduced. 


5.1.1. Lagrangian technigue 

By the Lagrangian moving grid method the governing equations (1) and (4) are 
solved by virtue of (34)-(37). The initial depth, of the element is taken to 
be the cell average of the exact initial profile. The initial velocity of the boundary, 
is given by the volume weighted velocity of the adjacent cells. They are 


= 

"'3 


IbO , 

J-1 


h(x, 0) dx 




h{x, 0) u{x, 0) dx 


hO _ hO 

°3 4-1 






/i(x, 0) dx 


(81) 


where bj and are the boundaries of the cell at t = 0, and c° denotes the 
initial centre of the cell. 

The constant inflow and outflow boundary conditions are executed by setting 
the depth gradient dh/dx at x = 5o(t) and x = bN{t) equal to zero, so that 
du/dt = 0 uo{t) = uo(0) and UN{t) = MAr(O) for < > 0 because the flow is on an 
non-accelerating slope Sx = 0. 

Fig. 7 demonstrates the simulated results {N = 60); oscillations develop as the 
shock wave passes through, and these persist even if the time step is selected to 
be very small. The velocities of the cell boundary after the shock are sometimes 
faster or slower than they should be and therefore oscillations take place. These 
oscillations propagate downslope as time increases and no shock wave propagates 
upslope. This indicates that the Lagrangian moving grid technique is ill behaved 
and cannot describe the travelling shock wave. 
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Lagrange, t=(p=6=40°, At=10 ^ Lagrange, t=(p=6=40°, At=10 ^ 



FIG. 7. Depth (left) and the corresponding velocity (right) profiles of the upslope travelling 
wave at t = 0, 3, 6, where circles denote the computed results at the cell centres and the solid line 
indicates the exact solution. The time step is taken to be At = 10“^ dimensionless time unit. 


NOCS (S), <^=0 = 6=40°, cfl = 0.3, Ax = 0.1 



FIG. 8. Depth profiles of the upslope travelling wave computed by the NOC scheme at 
t = 6 with N = 360. The solid lines indicate the exact solution and circles mean the computed 
results. 


5.1.2. Eulerian shock-capturing methods 

The NOC scheme is applied to (1) and (2) on a ID grid with 90 and 360 grid- 
points, respectively. The initial conditions are transferred to the mean values over 
the cells before the computing commences, 



'j-l/2 


■ 3 - 1/2 


(82) 








24 


Y.C. TAI, S. NOELLE, J.M.N.T. GRAY & K. HUTTER 


The constant inflow boundary condition is implemented by the assignments ho(t) = 
ho{0) and mo{t) = mo(0) at a; = 0. The outflow boundary condition is described 
by setting dhjdx = 0 and dm/dx = 0 a.t x = 36, where they are 

Un = (4 Un -1 - Un- 2 ) /3, for U = h, m, (83) 

by using the cell averages of the closest cells for a second order extrapolation. 

Three different cell reconstructions were tested: the NOC scheme with superbee 
limiter (NOCS-S), piecewise linear (r=2) and quadratic (r=3) WENO reconstruc¬ 
tions [26]. Fig. 8 demonstrates the simulated avalanche depth of the travelling wave 
problem (circles) and a comparison with the exact solution (solid line) at f = 6 di¬ 
mensionless time units. All of them are able to adequately describe this travelling 
shock wave problem. 


5.2. Parabolic similarity solution 

This section is concerned with the simulation of the parabolic similarity solution 
outlined in § 2.2. In the test problem the parabolic avalanche body is considered to 
slide on an inclined flat plane in the domain 0 < a; < 36 dimensionless length units 
with constant inclination angle C, = 40°. The basal and internal friction angles are 
simultaneously selected to be 30°, and the initial condition is chosen to be 50 = 1 
and po = 0. On the inclined plane the initial depth and velocity distributions are 
mapped into 


h{x, 0) = 1- ((x-4)/3.2)^ 
u(x, 0) = 1.2, 


for X G [0.8, 7.2]. 


(84) 


Our choice of the initial velocity, u(x, 0) = uq = 1.2, guarantees that condition 
(32) will be satisfied for all times. This problem will serve as the standard test 
problem for the resolution of the depth profile and the determination of the margin 
locations. 


5.3.1. Lagrangian technigue 

In the Lagrangian moving grid technique the model equations (1) and (4) are 
solved by virtue of (34)-(37) on a ID grid. The boundary condition is given 
by setting the heights at the margin (front and rear) points to be equal to zero, 
/io(x, t) = 0 and h^ix, t) = 0. 

Fig. 9 illustrates the simulated result at the dimensionless time units t = 0, 2, 4, 6 
with cell number N = 16, in which the circles denote the computed results and the 
solid line indicates the exact solution. The avalanche body extends as it flows 
down and still keeps the parabolic depth profile. The velocity is keeping a linear 
distribution through the bulk body. It ensures the symmetric depth profile during 
the motion. 

From the simulated results it follows that the Lagrangian moving grid technique 
can not only describe the depth profile well but also determines the margin locations 
of the similarity solution very accurately. There is excellent agreement between the 
simulated results and the exact solutions, see Fig. 9. The motions of the front and 
rear edges of the avalanche body in the similarity solution are illustrated in Fig. 10. 
The circles denote the computed results by the Lagrangian moving grid technique 
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6=(i,=30‘'. ^=40', At = 10-®. N=16 5=9S=30°, f=40", At=10-®. N=16 
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FIG. 9. Depth (left) and the corresponding velocity (right) profiles of the parabolic similarity 
solution (Problem I) computed by the Lagrangian moving-grid scheme at the dimensionless time 
units t = 0, 2, 4, 6, where the avalanche body is divided into 16 cells, and the time interval is 
At = 10-3. 





dimensionless distance, X 

FIG. 10. Locations of the front and rear edges of the avalanche body in the parabolic 
similarity solution problem as they evolve in time. The circles denote the computed results by 
the Lagrangian moving grid technique (A^ = 16), and the solid lines indicate the exact margin 
positions. They are in excellent agreement. 


and the solid lines indicate the exact locations of the margins. They are also in 
excellent agreement. 

The Lagrangian method is also tested by different grid numbers. Fig. 11 shows 
the results computed by different grid numbers, N = 16, 32, and 64, respectively. 
With different grid numbers this method can always keep the excellent resolutions 
when compared with the exact solutions. 

Caculations were also performed with initial condition pq ^ 0; results turned out 
to be similarly convincing as the above ones. For this reason they are not presented 
here [36]. 


5.2.2. Eulerian technique 
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FIG. 11. Depth profiles computed by the Lagrangian moving-grid technique for the parabolic 
similarity solution problem (Problem I), where the avalanche body is divided into different numbers 
of cells, N = 16, 32, 64. All the results are shown at t = 6 and the computational time interval is 
At = 10“^. The number of the cells does not influence the good agreement between the simulated 
results (circles) and the exact solutions (solid line). 


In § 3.2, the Eulerian schemes are based on the model equations (1) and (2) in 
conservative form, so that the velocity outside the avalanche body (inclusive the 
margin point) is not defined. Intuitively, adding a thin layer of the material over the 
whole computational domain could be used to treat the grain free regions. Another 
trick can also be introduced, in which all the physical variables are set to zero if 
h = 0. This would be reasonable since h = 0 —>■ m = hu = 0. 

Fig. 12 illustrates the comparison between the computed results obtained from 
the NOC scheme, where a thin layer Hq = lO”"* repectively ho = 0 is added over the 
whole computational domain, and from the scheme with our front-tracking method. 
All the three results of the depth profiles are acceptable except for the oscillation 
near the top. However, have a look at the velocity profiles in these figures, there 
are several cells with du/dx < 0 around the margins. This violates the assumption 
du/dx > 0 in the parabolic similarity solution problem. Moreover, the results show 
that there is large numerical diffusion around the margins (i.e. the margins move 
further than they should do) without the front-tracking method. For both reasons, 
the front-tracking method is needed to determine the location of the margins. 

Let us discuss the origin of the oscillation near the center of the avalanche. 
When one recomputes the solution using unlimited central differences for w', the 
oscillation disappears. Therefore, we have the following paradoxical situation: the 
introduction of TVD-limiters, which are needed to stabilize the solution in the 
presence of discontinuities, may destabilize the solutions in smooth regions! In fact, 
this is not entirely surprising, since in the presence of limiters the fluxes depend 
only Lipschitz-continuously on the data. 
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FIG. 12. Depth (left) and velocity (right) profiles of the parabolic similarity solution 

computed by the NOC scheme with Superbee limiter. In the top panels, a thin layer with ho = 
10“"^ is added to the whole computational domain. In the middle panels, all physical variables 
are set to zero if h = 0. Whilst, the bottom pannels demonstrate the results from the scheme with 
front-tracking method. The whole computational domain is divided into 90 cells (A^ = 90), the 
circles denote simulated results and the solid lines represent the exact solution. The results show 
that the added thin layer does not influence the depth profile very much, if it is sufficiently small, 
but the margin locations can not be exactly determined without the front-tracking method. An 
oscillation near the middle of the avalanche (local maximum) is visible in all three calculations. 
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FIG. 13. Depth (left) and velocity (right) profiles of the parabolic similarity solution 

at f = 6 computed by the NOC scheme with front-tracking and piecewise quadratic WENO 
cell reconstruction. The whole computational domain is divided into 90 cells {N = 90) and the 
Courant number is selected to be 0.3. The margin locations are well described and the oscillation 
near the center is successfully removed. 


We have therefore experimented with more smooth reconstructions, namely piece- 
wise quadratic WENO interpolants of Jiang and Shu [17] and Levy, Puppo and 
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FIG. 14. Front and rear edges of the avalanche body in the parabolic similarity solution 
simulated by the NOC-front-tracking scheme as they evolve in time, “o” denotes the computed 
results obtained with the piecewise quadratic WENO cell reconstruction, “x” means the results 
deduced with Superbee limiter and solid lines indicate the exact margin solution. 


Russo [26], which depend smoothly on the data and are at the same time nonoscil- 
latory at discontinuities. In the margin cells, we still use the piecewise linear re¬ 
constructions introduced in Section 4.1, and in the two cells adjacent to the margin 
cells, we use a piecewise linear WENO reconstruction. We have experimented both 
with second- and third-order quadrature rules in time. In our experience, both 
yield comparable results. Fig. 13 demonstrates the results for these reconstructions 
combined with our front-tracking method. The margin locations are well described 
by the front-tracking method, and the oscillation near the center is successfully 
removed (compare the bottom pannels in Fig. 12 and Fig. 13). 

Fig. 14 shows the computed front and rear edges of the avalanche body in the 
parabolic similarity solution as they evolve in time, “o” denotes the computed 
results obtained by the NOC scheme with the piecewise quadratic WENO cell 
reconstruction, “x” means the results deduced with Superbee limiter and solid lines 
indicate the exact margin solution. Both the Superbee limiter and the piecewise 
quadratic WENO cell reconstruction for the NOC front-tracking schemes can yield 
good agreement of the determined margin locations with the exact solutions. 

The use of the Superbee limiter results in a small delay of the avalanche body, i.e. 
a slower velocity both at the front and the rear. The reason is that the Superbee 
limiter tends to be overcompressive in smooth regions of the solution, and therefore 
it does not give the appropriate flux at the boundaries between the internal and 
the margin cells. 

In order to obtain some quantitative information on the accuracy of the schemes, 
we introduce an error measure for the depth by 


N 

E 

j=0 


-j-exact 
hj hj 


—exact 

j=o 


E = 


(85) 
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TABLE 5.1 

Error (85) of the different schemes. 



NFT(S90) 

NFT(W90) 

NFT(S360) 

NFT(W360) 

Lag(16) 

Lag(32) 


E(xl0"®) 

E(xl0“®) 

t = 

1 

9.0247 

7.6096 

0.8816 

0.8813 

1.7130 

0.2937 

t = 

2 

11.3333 

7.5679 

0.8532 

0.9024 

1.7764 

0.3664 

t = 

3 

12.8854 

7.1051 

0.7336 

0.9492 

1.8944 

0.4135 

t = 

4 

16.3450 

6.0860 

0.8484 

1.0875 

1.8888 

0.4413 

t = 

5 

17.4274 

6.5434 

0.8672 

1.1203 

1.8974 

0.4492 

t = 

6 

18.8426 

6.5340 

1.1109 

1.3203 

1.9474 

0.4658 

t = 

7 

16.5043 

6.0608 

1.3970 

1.0435 

1.8817 

0.5026 

t = 

8 

15.3966 

6.0762 

1.8558 

1.1981 

1.9526 

0.4830 


where denotes the cell averaged depth of the exact solution. The errors 

of the Lagrangian method, the NOC front-tracking scheme with Superbee limiter 
(NFT(S)) and piecewise quadratic WENO interpolations (NFT(W)) at t = 1 to 8 
dimensionless time units are shown in Table. Here, the Eulerian schemes are tested 
by using N = 90 and N = 360, respectively, and for the Lagrangian scheme N = 16 
and N = 32 are used. The Lagrangian method results in the least errors, obviously 
smaller than all the Eulerian schemes. It also converges at a better rate. 

5.3. Upward Moving Shock Wave 

Shock formations are often observed when the avalanche slides into the run-out 
horizontal zone. Here the front part comes to rest, while the tail accelerates further 
and its velocity becomes supercritical. In [38] a comparison was made between our 
shock-capturing method and the Lagrangian moving grid technique for the case of 
coinciding basal and internal friction angles. Here we compute a flow with basal 
friction angle (j) = 38° and internal friction angle <5 = 35°. As a consequence, we 
have a jump in the earth-pressure coefficient when the flow changes from an 
expanding {ux > 0) to a contracting region (ux < 0). 

The setup is as follows: The granular material released from a parabolic cap 
slides down an inclined plane and merges into the run-out horizontal zone. The 
centre of the cap is initially located at a; = 4.0 and the initial radius and the height 
are 3.2 and 1.0 dimensionless length units, respectively. The inclination angle of the 
inclined plane is 40° and the (linear and continuous) transition region lies between 
X = 21.5 and x = 25.5. We use 180 gridpoints and a CFL-number of 0.4. 

Fig. 15 illustrates the simulated process as the avalanche slides on the inclined 
plane into the horizontal run-out zone (so initially the flow is expanding). The 
avalanche body extends on the inclined plane until the front reaches the run-out 
zone. Here the basal friction is enough to bring the front of the granular material to 
rest while the rear part accelerates further. Therefore, the flow becomes contracting 
in the transition zone. At this stage, a shock (surge) wave is created {t = 12), 
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FIG. 15. Process of the avalanche simulated by the shock-capturing and front-tracking 

NOC method at t = 0, 3, 6,..., 27 dimensionless time units. As the front reaches the run-out 
zone and comes to rest, the part of the tail accelerates further and the avalanche body contracts. 
Once the velocity becomes supercritical, a shock wave develops, which moves upward. The dashes 
below the graphs mark the tail and the head of the avalanche. 


which moves upward. Such shock waves make the Lagrangian method unstable, if 
no artificial viscosity is applied (see [38]). Our non-oscillatory central front-tracking 
scheme handles both the shock wave and the margins of the avalanche well. 

6. CONCLUSION 

In this paper we have developed a Lagrangian and an Eulerian shock-capturing 
finite-difference scheme with front-tracking for the spatially one-dimensional Savage- 
Hutter equations of granular avalanches. The purpose was to reproduce the tem¬ 
poral evolution of the avalanche geometry and downslope velocity under situations 
when internal shocks may occur. This happens e.g. when an avalanche of finite 
mass moves from an inclined chute into the horizontal run-out zone and, in the 
transition zone, is deccelerated from a supercritical flow state to a subcritical state. 
The Lagrangian scheme (which is excellent for smooth solutions) develops unphys¬ 
ical oscillations when the solution contains, or develops, shock discontinuities. In 
order to compute discontinuous solutions, we propose to use a conservative shock¬ 
capturing finite difference scheme. We adapt the second-order accurate staggered 
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scheme of Nessyahu and Tadmor [30] to the Savage-Hutter equations. The stag¬ 
gered approach avoids the use of characteristic decompositions which are needed in 
standard upwind schemes, but are not known for the Savage-Hutter equations. We 
show that our non-oscillatory central (NOC) scheme reproduces both smooth and 
shock solutions adequately except for the following two problems: First, oscillations 
may occur near smooth extrema due to the presence of piecewise linear reconstruc¬ 
tions with TVD-type limiters. These oscillations disappear when one uses piecewise 
quadratic cell reconstructions in the interior of the avalanche. Second, our NOC 
scheme (and in fact, any Eulerian scheme) does not capture the vacuum-boundary 
accurately. This may lead to serious stability problems. We improve the treat¬ 
ment of the free boundary by combining the scheme with a front-tracking method 
applied to the margin cells. In the spirit of the Nessyahu-Tadmor scheme, we do 
not make use of the vacuum Riemann-problem, but rely on a new piecewise linear 
reconstruction at the vacuum boundary and carefully chosen Taylor-extrapolations 
for the corresponding numerical fluxes. With such a combination of an internal 
Eulerian NOC scheme and a Lagrangian “boundary scheme” two standard test 
problems - an upward moving shock and a parabolic cap moving down an inclined 
plane - could be well reproduced. The scheme also produces satisfactory results for 
the more realistic problem mentioned above: an avalanche moving down an inclined 
plane and coming to rest at a flat run-out. Here an upward moving shock wave 
develops from smooth data, and the flow changes from being expanding to con¬ 
tracting ahead of the shock. In this situation, the earth pressure coefficient changes 
discontinuously, so we are facing the full difficulties inherent in the Savage-Hutter 
model. 

Several questions remain and await further study: 

• The shock-capturing NOC numerical method including the front-tracking scheme 
should be extended to two-dimensional flows. This is work under progress. 

• The original Lagrangian moving grid scheme could also be developed as a 
shock-capturing scheme. Here the main difficulty would be in the determination of 
the correct grid velocity. 

We are working on these topics and will report on results in due time. 
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